06KASP流程-通用-有参重测序&简化 or 无参简化
一.有参kasp
位点交付文件
kasp_analysis/conventional_primer/samples.final.KASP.info.xls
指纹图等重测序或简化流程生成vcf文件再做
# 以analysis/vcf_filter/samples.pop.snp.recode.vcf为例
# 示例是analysis/variation_dir/variants_anno_dir 太大
# analysis/vcf_filter/samples.gatk.snp.m2M2.mis0.2.mac3.4dtv.vcf 备用
1.1 配置文件
/share/nas2/pub/pipline/dna-seq/kasp/v1.0/kasp-develop.config.yaml
mkdir kasp
cd kasp
cp /share/nas6/zhouxy/project/ddrad/GP-20220524-4388-1_guizhoudaxue_chentaolin_48samples_ddrad/kasp-develop.config.yaml kasp-develop.config.yaml
修改配置文件
1.2 主流程
/share/nas6/zhouxy/pipline/kasp/kasp-develop/v1.2/kasp-develop.pl
/share/nas1/yuj/pipline/kasp/kasp-develop/v1.2/kasp-develop.pl
/share/nas2/pub/pipline/dna-seq/kasp/v1.0/kasp-develop.pl
/share/nas2/pub/pipline/dna-seq/kasp/v1.0/kasp-develop.v1.1.pl # 最新
perl /share/nas1/yuj/pipline/kasp/kasp-develop/v1.2/kasp-develop.pl -i kasp-develop.config.yaml -o kasp_analysis
perl /share/nas2/pub/pipline/dna-seq/kasp/v1.0/kasp-develop.v1.1.pl -i kasp-develop.config.yaml -o kasp_analysis # 先用这个试试
下面这个脚本是等待vcf生成就开始kasp流程
#!/bin/bash
check_file_exists() {
if [[ -f "$1" ]]; then
return 0 # 文件存在
else
return 1 # 文件不存在
fi
}
# 等待文件生成
while ! check_file_exists "/share/nas1/yuj/project/re-sequencing/GP-20230517-6067_20230621/analysis_dir/variation_dir/variants_anno_dir/samples.pop.snp.anno.result.vcf"; do
sleep 60 # 每秒检查一次文件是否存在
done
echo "ann.xls已生成，执行A命令"
nohup perl /share/nas6/zhouxy/pipline/kasp/kasp-develop/v1.1/kasp-develop.pl -i kasp-develop.config.yaml -o kasp_analysis &
1.3 整理结果
## 原始v1.1路径
# perl /share/nas6/zhouxy/pipline/kasp/kasp-develop/v1.1/bin/resultsdir/resultsdir_v2.pl -i kasp_analysis/ -o kasp_complete_dir
perl /share/nas1/yuj/pipline/kasp/kasp-develop/v1.2/bin/resultsdir/resultsdir_v2.pl -i kasp_analysis/ -o kasp_complete_dir
1.4 生成报告
cp /share/nas6/zhouxy/pipline/kasp/kasp-develop/v1.1/html_report/kasp_report.cfg . && realpath kasp_report.cfg
修改配置
## 原路径
# perl /share/nas6/zhouxy/pipline/kasp/kasp-develop/v1.1/html_report/kasp_Web_Report.pl -id kasp_complete_dir/ -cfg kasp_report.cfg
perl /share/nas1/yuj/pipline/kasp/kasp-develop/v1.2/html_report/kasp_Web_Report.pl -id kasp_complete_dir/ -cfg kasp_report.cfg
二.FAQ
2.1 注意事项及问答
2.1.1 注意事项
1.关于KASP的标记选择，要把观测杂合度设置为0.3以下，超过0.3的标记不要给。如果是0-0.2就更好了
2.在我们所有GATK callsnp结果中，都会有一种情况，就是明明AD等位基因深度为0,0，但是基因型是0/0，此种情况是软件导致的。
我们的应对方法是：把深度低于3X（包括1X,2X，但是不包括3X）的等位基因型归为./.（位置），此问题使用vcftools过滤处理，加参数--minDP 3。或者自己写程序修改。
vcftools在过滤后生成的结果文件中，把小于深度3的输出为./.（参数--minDP 3）。
所以我们原始结果的一些统计是有问题的。需要修复下。没有深度覆盖的基因型是有问题的。
目前涉及的流程有重测序，简化测序，群体进化等使用GATK软件的流程，对于BSA,BSR影响较小，BSA是对每个样品单独过滤的。
2.1.2 结果中03/04/05文件夹内标记变化的原因
首先啊，客户看的是第一次的结果, 第一次的结果不能用啊。
-----------------------------------------------------------
我拿第二次的结果来说啊，第二次是剔除了污染的样品：
03文件夹里是引物设计成功的201个，然后这201个与一代引物设计的结果取了交集，
就变成了04文件夹里的77个。
然后对这77个，每一个位点都计算了其对样品的区分能力，区分能力从大到小排了个序，
然后截取前50个放在05文件夹里，与此同时也计算了区分全部样品最少只需要10个就够了。
50个是77个排序之后的前50个
这10个是77个排序之后的前10个，也是这50个的前10个
------------------------------------------------------------
再说回第一次结果也是这样的，不过是03文件夹94个，04文件夹里就40个了
那到05里头取前50个也只能有40个
2.1.3 一代引物设计怎么来的
"前面开发引物成功"指的是 只要在snp位点两侧能设计成功引物就行
这个所谓的"一代引物"是在snp位点邻近的两侧还有ssr标记,且设计引物能把包括snp位点和ssr标记都跑出来
<h5 id="a42" class="title-h5">4.2.3 SSR标记辅助筛选</h5>
<p class="paragraph">在成功设计的KASP标记基础上，进一步在SNP位点上下游500bp范围内系统搜索SSR（Simple Sequence Repeat）重复序列。筛选标准要求：</p>
<div class="paragraph">
<p>1. SSR重复单元类型要求：</p>
<ul>
<li>单核苷酸重复：≥4次重复</li>
<li>二核苷酸重复：≥4次重复</li>
<li>三核苷酸重复：≥6次重复</li>
<li>四核苷酸重复：≥5次重复</li>
<li>五核苷酸重复：≥4次重复</li>
<li>六核苷酸重复：≥4次重复</li>
</ul>
<p>2. SSR重复区域需满足以下条件：</p>
<ul>
<li>距离SNP位点>50bp</li>
</ul>
</div>
<p class="paragraph">对符合条件的SSR标记进行PCR引物设计，最终筛选出1个可同时扩增SNP位点和SSR标记的复合引物，具体设计结果请查看文件夹04_selectMarker。</p>
<li><a href="#a33" class="js-active-bg">4.2 SNP位点过滤和引物设计</a>
<ul class="nav">
<li><a href="#a34" class="js-active-bg">4.2.1 SNP位点步骤</a></li>
<li><a href="#a41" class="js-active-bg">4.2.2 KASP位点引物设计</a></li>
<li><a href="#a42" class="js-active-bg">4.2.3 SSR标记辅助筛选</a></li>
</ul>
</li>
2.2 叶绿体项目
2.2.1 vcf_filter.cmds没有位点
其中第二个程序重新运行
1）完全解除过滤
--max-missing 1.0 \ # 允许100%缺失
--min-alleles 1 \ # 允许单等位
--max-alleles 100 \ # 允许多等位
--maf 0.0 \ # 关闭MAF过滤
--minQ 0 \ # 关闭质量过滤
--min-meanDP 0 \ # 关闭深度过滤
2）ai给出的建议
--min-alleles 1 \ # 允许单等位
--max-alleles 1 \ # 限制为单等位
--min-meanDP 1 \ # 最低平均深度更宽松
--maxDP 10000 \ # 适当提高最大深度阈值
--minQ 20 \ # 降低质量值阈值
--max-missing 0.8 \ # 允许更高缺失率
--maf 0.0 # 关闭MAF过滤
3）成功运行的参数
/share/nas6/zhouxy/biosoft/vcftools/current/bin/vcftools --vcf kasp/kasp_analysis/vcf_filter/samples.snp.region.filter.vcf --recode --out kasp/kasp_analysis/vcf_filter/samples.snp.filter --min-alleles 1 --max-alleles 2 --min-meanDP 0 --maxDP 10000 --minQ 20 --max-missing 0.9 --maf 0.0 --minDP 3
2.2.2 PIC_filter.cmds
perl /share/nas6/zhouxy/pipline/kasp/kasp-develop/v1.2/bin/stats/kaspfilter_Marker2stat_line.pl -i /share/nas2/yuj/project/2025/chloroplast/GP-20211022-3553-6_20250210/data/kasp/kasp_analysis/blast_filter/samples.snp.blast.filter.xls -o /share/nas2/yuj/project/2025/chloroplast/GP-20211022-3553-6_20250210/data/kasp/kasp_analysis/PIC_filter -nrtio 0 -ho 0.6 -he 0.2 -pic 0.30
脚本会保留满足以下条件的SNP：
观测杂合率（实际杂合样本比例）≤ 0.6
期望杂合率（基于等位基因频率的理论值）≥ 0.2
PIC值（多态性信息含量）≥ 0.30
非参考型等位基因比例未被限制（-nrtio 0 表示关闭此过滤）
perl /share/nas6/zhouxy/pipline/kasp/kasp-develop/v1.2/bin/seqalign/region_seq2substr_merge.pl -r /share/nas2/yuj/project/2025/chloroplast/GP-20211022-3553-6_20250210/data/kasp/kasp_analysis/PIC_filter/samples.snp.PIC.filter.stat.xls -g /share/nas2/yuj/project/2025/chloroplast/GP-20211022-3553-6_20250210/data/kasp/kasp_analysis/genome/genome.fa -o /share/nas2/yuj/project/2025/chloroplast/GP-20211022-3553-6_20250210/data/kasp/kasp_analysis/PIC_filter/samples.snp.PIC.seq.xls
2.3 vcf_filter结果太少
可能的原因
1.因为没有注释，无法筛选cds上的，选择扩大间隔，-l 500
1.
perl /share/nas6/zhouxy/pipline/kasp/kasp-develop/v1.2/bin/data_access/get_50bp_list.pl -i analysis_dir/vcf_filter/samples.pop.snp.filter.format.recode.vcf -o kasp/kasp_analysis/vcf_filter/samples.snp.region.filter.vcf -l 500
# 间隔默认-l 100
2.
## 常规2倍体项目
/share/nas6/zhouxy/biosoft/vcftools/current/bin/vcftools --min-meanDP 5 --maxDP 500 --minDP 3 --minQ 30 --max-missing 0.90 --maf 0.05 --min-alleles 2 --max-alleles 2 --vcf kasp/kasp_analysis/vcf_filter/samples.snp.region.filter.vcf --recode --out kasp/kasp_analysis/vcf_filter/samples.snp.filter
# --minDP 3 滤去没有覆盖度的位点
3.
perl /share/nas6/zhouxy/pipline/kasp/kasp-develop/v1.2/bin/snp2xxx/vcf_to_snplist_v1.2.pl -i kasp/kasp_analysis/vcf_filter/samples.snp.filter.recode.vcf -o kasp/kasp_analysis/vcf_filter/samples.snp.filter.recode.tab.xls -ref 1
2.4 指纹图显示不全
1.第二种样式
perl /share/nas6/xul/program/ddrad/fingerprint_info_heatmap.pl -i samples.fingerprint.info.xls # 另一种样式
2.
修改/share/nas6/zhouxy/pipline/kasp/kasp-develop/v1.1/bin/fingerprint/plot_pheatmap.r
复制到自己目录一份/share/nas1/yuj/pipline/kasp/kasp-develop/v1.1/bin/fingerprint/plot_pheatmap.r
修改其中的pdf(args[2],50,50)后两个参数表示宽度和高度
对应kasp_analysis/.cmds_dir/selectVariant.cmds的最后两个命令
Rscript /share/nas1/yuj/pipline/kasp/kasp-develop/v1.1/bin/fingerprint/plot_pheatmap.r kasp_analysis/fingerprint/samples.fingerprint.info.plot.xls kasp_analysis/fingerprint/samples.fingerprint.info.plot.pdf && convert -density 600 kasp_analysis/fingerprint/samples.fingerprint.info.plot.pdf kasp_analysis/fingerprint/samples.fingerprint.info.plot.png
2.5 extract.cmds中gff3ToGenePred报错
2.5.1 类型1
1.报错信息
Annotation records for discontinuous features with ID="TRNAQ-CUG-7" do not have the same type, found "tRNA" and "gene"
Annotation records for discontinuous features with ID="TRNAQ-CUG-8" do not have the same type, found "tRNA" and "gene"
Annotation records for discontinuous features with ID="TRNAQ-CUG-9" do not have the same type, found "tRNA" and "gene"
Annotation records for discontinuous features with ID="TRNAD-GUC-10" do not have the same type, found "tRNA" and "gene"
GFF3: 51 parser errors
2.显示一下gff第三列
cp kasp_analysis/genome/genome.gff3 kasp_analysis/genome/genome.gff3.bak
cut -d
### 2.5.2 类型2
```shell
1.报错信息
Annotation records for discontinuous features with ID="Capann_59Chr02g031500.1.five_prime_UTR" do not have the same type, found "three_prime_UTR" and "five_prime_UTR"
Annotation records for discontinuous features with ID="Capann_59Chr02g031500.1.five_prime_UTR" do not have the same type, found "three_prime_UTR" and "five_prime_UTR"
Annotation records for discontinuous features with ID="Capann_59Chr02g031490.1.five_prime_UTR" do not have the same type, found "three_prime_UTR" and "five_prime_UTR"
GFF3: 51 parser errors
id中的类型与实际类型不一样 前后不一致
2.显示一下gff第三列
cp kasp_analysis/genome/genome.gff3 kasp_analysis/genome/genome.gff3.bak
cut -d
### 2.5.3 类型3 无效的属性标签
```shell
invalid attribute tag, must start with an alphabetic character and be composed of alphanumeric or underscore characters: dev-stage
确保属性标签以字母开头，并只包含字母、数字或下划线。
就把gff3中dev-stage全部替换成dev_stage就行
2.6 基因型文件染色体名称替换
awk 'BEGIN {OFS="\t"}
# 处理映射文件（使用空格分隔）
NR==FNR {chrom_map[$1]=$2; next}
# 处理数据文件（切换为制表符分隔）
FNR==1 {FS=OFS; print; next}
{
if ($1 in chrom_map) {
$1 = chrom_map[$1]
}
print
}' chromosome_mapping.txt samples.kasp.genotype.select.xls > samples.kasp.genotype.select.renamed.xls
2.7 挑选指定区间基因型
python3 /share/nas1/yuj/script/ddradANDreseqANDgenome/introgressed_fragment/find_overlapping_intervals.py --files /share/nas2/yuj/project/2025/re-sequencing/GP-20241021-9521-2_20250104/02daerwenshimian/20250325plot/fig1/125K_vcf_stat_vs_refgenome/group1/*.txt # 生成共有区间
python3 /share/nas1/yuj/script/ddradANDreseqANDgenome/introgressed_fragment/interval_genotype_filter.py --input samples.kasp.genotype.select.renamed.xls --intervals overlapping_intervals.txt --output filtered_results.xls # 挑出共有区间的基因型数据
三.回推鉴定
cd xx
控制selectVariant/samples.kasp.genotype.select.xls里面的对应位点
perl /share/nas6/zhouxy/pipline/kasp/kasp-develop/v1.2/bin/fingerprint/findMinimumMarker.pl -i kasp_analysis/selectVariant/samples.kasp.genotype.select.xls -o kasp_analysis/fingerprint -n 50 && perl /share/nas6/zhouxy/pipline/kasp/kasp-develop/v1.2/bin/fingerprint/markerIdent.pl -i kasp_analysis/fingerprint/samples.sortMarker.xls -o kasp_analysis/fingerprint/samples.sortMarker.ident.xls && Rscript /share/nas6/zhouxy/pipline/kasp/kasp-develop/v1.2/bin/fingerprint/ideCoeff.r kasp_analysis/fingerprint/samples.sortMarker.ident.xls kasp_analysis/fingerprint/samples.sortMarker.ident.pdf `perl -lane 'next if(/Num/);if($F[1]>=1){print $F[0];exit;}' kasp_analysis/fingerprint/samples.sortMarker.ident.xls` && convert -density 600 kasp_analysis/fingerprint/samples.sortMarker.ident.pdf kasp_analysis/fingerprint/samples.sortMarker.ident.png
四.无参简化指纹图
主流程：
perl /share/nas6/zhouxy/pipline/kasp/kasp-develop/v1.1/kasp-develop_ddrad_with_noref.pl -v samples.pop.snp.recode.vcf -g ../stacks_stat/consensus/tags.consensus.fa
1，按照深度(10x）、标记完整度（0。9）、maf（0.05）过滤
2，按照比对NR库上的序列过滤，保留比对上NR库的
3，按照PIC（0.35）过滤。
4，kasp引物设计
5，指纹图
注：如果第一步后标记很多，第二步后标记很少，可以增加第一步参数，跳过第二步。
整理结果：
perl /share/nas6/zhouxy/pipline/kasp/kasp-develop/v1.1/bin/resultsdir/resultsdir_ddrad_with_noref.pl -i kasp_analysis -o kasp_result
报告：
cp /share/nas6/zhouxy/pipline/kasp/kasp-develop/v1.1/html_report/kasp_report.cfg . && realpath kasp_report.cfg
perl /share/nas6/zhouxy/pipline/kasp/kasp-develop/v1.1/html_report/kasp_ddrad_with_noref_Web_Report.pl -id kasp_result/ -cfg kasp_report.cfg
五.流程详解
# 初始化阶段
##----------------------------------------------------------------------------------------------
# 读取配置文件，设置工具路径（samtools/blastn/primer3等）
my %config=&readConfig("$Bin/bin/Config/CFG");
# Step 00 : 配置检查
# 读取物种、前缀、VCF、GFF、序列文件等基础配置
# Step 01 : 参考基因组处理
# 建立基因组索引，生成染色体长度文件，替换VCF中的序列名称
# Step 02 : 信息提取
# 转换GFF3为GenePred格式，提取mRNA序列用于后续注释
# Step 03 : 变异注释
# 使用ANNOVAR进行变异注释，生成多维度注释结果
# Step 04 : VCF过滤
# 通过质量值、缺失率、MAF等参数过滤原始SNP
# 生成过滤后的VCF和SNP列表
# Step 05 : 序列比对过滤
# 提取SNP侧翼序列进行blast比对，排除多比对位点
# Step 06 : PIC值过滤
# 基于多态信息含量(PIC)筛选高质量标记
# 进行群体遗传分析（杂合度、遗传多样性等）
# Step 07 : 引物设计
# 使用primer3进行KASP引物设计
# 通过两轮blast筛选特异性引物：
# 1. 初步过滤多比对引物
# 2. 严格筛选唯一比对引物
# 设计常规验证引物并进行类似过滤
# Step 10 : 指纹图谱分析
# 筛选最终标记组合
# 进行个体识别率分析和热图可视化
# Step 11 : 系统发育树
# 基于SNP数据构建最大似然树
# Step 12 : 群体结构分析
# 使用structure软件进行群体遗传结构分析
# Step 13 : PCA分析
# 主成分分析展示样本间遗传关系
# 辅助函数
##----------------------------------------------------------------------------------------------
sub Pipline_sh_commands { # 任务调度核心方法
# 使用Pipeliner模块管理任务依赖关系
# 通过flag文件跟踪任务状态
}
sub readConfig { # 配置文件读取方法
# 解析CFG配置文件获取工具路径
}
\t' -f 3 kasp_analysis/genome/genome.gff3.bak | sort | uniq | grep -v "#"
cDNA_match
CDS
exon
gene
lnc_RNA
mRNA
pseudogene
region
rRNA
snoRNA
snRNA
transcript
tRNA
3.正常的gff
CDS
exon
gene
mRNA
region
4.解决办法xul
perl -lane 'if(/\tGnomon\t/){print}elsif(/#/){print}' /share/nas6/pub/genome/shizi/Diospyros_lotus/GCF_014633365.1/GCF_014633365.1_ASM1463336v1_genomic.gff > new.gff
### 2.5.2 类型2
{{CODE_BLOCK_15}}
### 2.5.3 类型3 无效的属性标签
{{CODE_BLOCK_16}}
## 2.6 基因型文件染色体名称替换
{{CODE_BLOCK_17}}
## 2.7 挑选指定区间基因型
{{CODE_BLOCK_18}}
# 三.回推鉴定
{{CODE_BLOCK_19}}
# 四.无参简化指纹图
{{CODE_BLOCK_20}}
# 五.流程详解
{{CODE_BLOCK_21}}
\t' -f 3 kasp_analysis/genome/genome.gff3.bak | sort | uniq | grep -v "#"
CDS
exon
gene
mRNA
region
five_prime_UTR
three_prime_UTR
3.正常的gff
CDS
exon
gene
mRNA
region
4.解决办法
sed 's/three_prime_UTR/five_prime_UTR/g' /share/nas1/yuj/project/ddRAD-seq/GP-20230814-6613_20231213/data/Ca_59.v0.gff3 > /share/nas1/yuj/project/ddRAD-seq/GP-20230814-6613_20231213/data/new.gff
然后修改配置文件中gff路径
重新运行
2.5.3 类型3 无效的属性标签
{{CODE_BLOCK_16}}
2.6 基因型文件染色体名称替换
{{CODE_BLOCK_17}}
2.7 挑选指定区间基因型
{{CODE_BLOCK_18}}
三.回推鉴定
{{CODE_BLOCK_19}}
四.无参简化指纹图
{{CODE_BLOCK_20}}
五.流程详解
{{CODE_BLOCK_21}}
\t' -f 3 kasp_analysis/genome/genome.gff3.bak | sort | uniq | grep -v "#"
cDNA_match
CDS
exon
gene
lnc_RNA
mRNA
pseudogene
region
rRNA
snoRNA
snRNA
transcript
tRNA
3.正常的gff
CDS
exon
gene
mRNA
region
4.解决办法xul
perl -lane 'if(/\tGnomon\t/){print}elsif(/#/){print}' /share/nas6/pub/genome/shizi/Diospyros_lotus/GCF_014633365.1/GCF_014633365.1_ASM1463336v1_genomic.gff > new.gff
### 2.5.2 类型2
{{CODE_BLOCK_15}}
### 2.5.3 类型3 无效的属性标签
{{CODE_BLOCK_16}}
## 2.6 基因型文件染色体名称替换
{{CODE_BLOCK_17}}
## 2.7 挑选指定区间基因型
{{CODE_BLOCK_18}}
# 三.回推鉴定
{{CODE_BLOCK_19}}
# 四.无参简化指纹图
{{CODE_BLOCK_20}}
# 五.流程详解
{{CODE_BLOCK_21}}